#include "calculate.h"
#include <cmath>
#include <iostream>

using std::complex;

//返回单个域的相对介电常数
double d_perm(Para para, int j)
{
    int index = (j-1) % para.get_diel_n();
    double perm[MAX];
    para.get_diel_perm(perm);
    return perm[index];
}
//得到所有域的相对介电常数取值数组
void d_perm(Para para, double _epsilon[], int n)
{
    _epsilon[0] =1.0;
    _epsilon[n+1] = 1.0;
    for (int i = 1; i <= n; i++)
    {
        _epsilon[i] = d_perm(para, i);
    }
}


//返回单个域内的kz取值
double kz(double _epsilon[], double _omega, double _theta, int j)
{
    return _omega / C *sqrt(_epsilon[j] - sin(_theta) * sin(_theta));
}
//得到所有域的kz的取值数组
void kz(double _epsilon[], double _omega, double _theta, double k[], int n)
{
    for (int i = 0; i < n+2; i++)
    {
        k[i] = kz(_epsilon, _omega, _theta, i);
    }
}


//计算单个界面的反射系数
double R(double k[], double _epsilon[], int dir, int j)
{
    double temp1 = k[j]/_epsilon[j]; double temp2 = k[j+1]/_epsilon[j+1];

    if (dir == 1)
    {
        return (temp1 - temp2)/(temp1 + temp2);
    } else {
        return -(temp1 - temp2)/(temp1 + temp2);
    }
}
//得到所有界面的前向反射系数
void R1(double k[], double _epsilon[], double Rf[], int n)
{
    for (int i = 0; i < n+1; i++)
    {
        Rf[i] = R(k, _epsilon, 1, i);
    }
}
//得到所有界面的后向反射系数
void R2(double k[], double _epsilon[], double Rb[], int n)
{
    for (int i = 0; i < n+1; i++)
    {
        Rb[i] = R(k, _epsilon, -1, i);
    }
}


//计算单个界面的透射系数
double T(double k[], double _epsilon[], int dir, int j)
{
    double kz1 = k[j];
    double kz2 = k[j+1];
    double _epsilon1 = _epsilon[j];
    double _epsilon2 = _epsilon[j+1];
    double temp1 = kz1/_epsilon1; double temp2 = kz2/_epsilon2;

    if (dir == 1)
    {
        return 2*temp1/(temp1+temp2) * _epsilon1/_epsilon2;
    } else {
        return 2*temp2/(temp1+temp2) * _epsilon2/_epsilon1;
    }
}
//得到所有界面的前向透射系数
void T1(double k[], double _epsilon[], double Tf[], int n)
{
    for (int i = 0; i < n+1; i++)
    {
        Tf[i] = T(k, _epsilon, 1, i);
    }
}
//得到所有界面的后向反射系数
void T2(double k[], double _epsilon[], double Tb[], int n)
{
    for (int i = 0; i < n+1; i++)
    {
        Tb[i] = T(k, _epsilon, -1, i);
    }
}


//计算所有界面的全局前向反射系数
void GR1(Para para, const double k[], const double Tf[], const double Rf[], const double Tb[], const double Rb[], complex<double> GRf[])
{
    int n = para.get_diel_n() * para.get_r_p_n();
    double l[para.get_diel_n()]; para.get_diel_l(l);
    int index;
    double d;

    GRf[n] = complex<double>(Rf[n]);

    for (int i = n-1; i >= 0; i--)
    {
        index = para.get_diel_n() - (n-i)%para.get_diel_n();
        d = l[index];
        complex<double> c(cos(2*k[i+1]*d*um),sin(2*k[i+1]*d*um));
        GRf[i] = Rf[i] + (Tf[i]*Tb[i]*GRf[i+1]*c)/(1.0-Rb[i]*GRf[i+1]*c);
    }
}
//计算所有界面的全局后向反射系数
void GR2(Para para, const double k[], const double Tf[], const double Rf[], const double Tb[], const double Rb[], complex<double> GRb[])
{
    int n = para.get_diel_n() * para.get_r_p_n();
    double l[para.get_diel_n()]; para.get_diel_l(l);
    int index;
    double d;

    GRb[0] = complex<double>(Rb[0]);

    for (int i = 1; i <= n; i++)
    {
        index = (i-1)%para.get_diel_n();
        d = l[index];
        complex<double> c(cos(2*k[i]*d*um),sin(2*k[i]*d*um));
        GRb[i] = Rb[i] + (Tf[i]*Tb[i]*GRb[i-1]*c)/(1.0-Rf[i]*GRb[i-1]*c);
    }
}


//计算因子N
void facN(Para para, const double k[], const double Tf[], const double Rf[], const complex<double> GRb[], complex<double> _N[])
{
    int n = para.get_diel_n() * para.get_r_p_n();
    double l[para.get_diel_n()]; para.get_diel_l(l);
    int index;
    double d;

    for (int i = 0; i < n; i++)
    {
        index = i%para.get_diel_n();
        d = l[index];
        complex<double> c(cos(k[i+1]*d*um), sin(k[i+1]*d*um));
        _N[i] = c * Tf[i+1] / (1.0 - Rf[i+1] * GRb[i] * c * c);
    }
}
//计算每个域j总辐射场中同N有关的因子项
void relatedN(int j, const int n, const complex<double> _N[], complex<double> _rN[])
{
      if(j > 2 && j <= n+2)
      {
          _rN[j-2] = 1.0;
          for(int k = j-3; k >=0; k--)
          {
              _rN[k] = _N[k] *_rN[k+1];
          }
      } else
      {
          std::cout << "因子N计算相关参数传递错误" << std::endl;
      }
}


//计算因子O
void facO(Para para, const double k[], const double Tb[], const double Rb[], const complex<double> GRf[], complex<double> _O[])
{
    int n = para.get_diel_n() * para.get_r_p_n();
    double l[para.get_diel_n()]; para.get_diel_l(l);
    int index;
    double d;

    for (int i = 0; i < n; i++)
    {
        index = i%para.get_diel_n();
        d = l[index];
        complex<double> c(cos(k[i+1]*d*um), sin(k[i+1]*d*um));
        _O[i] = c * Tb[i] / (1.0 - Rb[i] * GRf[i+1] * c * c);
    }
}
//计算域j总辐射场中与因子O相关的连乘数组
void relatedO(int j, const int n, const complex<double> _O[], complex<double> _rO[])
{
    if(j < n+1 && j >0)
    {
        _rO[0] = 1.0;
        for(int k = 1; k < n-j+2; k++)
        {
            _rO[k] = _rO[k-1]*_O[j-2+k];
        }

        for(int k = 0; k < n-j+1; k++)
        {
            _rO[k] = _rO[k+1];
        }
    } else 
    {
        std::cout << "因子O计算相关参数传递错误" << std::endl;
    }
}


//计算每层向前辐射的振幅
void source1(const double _omega, const double par_v, Para para, const double k[], const double _epsilon[], complex<double> sf[])
{
    int n = para.get_diel_n() * para.get_r_p_n();
    double l[para.get_diel_n()]; para.get_diel_l(l);
    int index;
    double d = 0.0;
    double temp = _omega * _omega / (C * C) - k[0] * k[0];
    double num1, num2; double den1, den2, den3;
    complex<double> c1, c2;

    for(int i = 0; i < n+1;i++)
    {
        num1 = par_v/C*(temp*C*C)/(_omega*_omega*_epsilon[i+1])*(_epsilon[i+1]-_epsilon[i]);
        num2 = 1-par_v*par_v/(C*C)*_epsilon[i+1]-par_v/C*k[i]/(_omega/C);
        den1 = 1-par_v*par_v/(C*C)*_epsilon[i+1]+temp*par_v*par_v/(_omega*_omega);
        den2 = 1-par_v/C*k[i]/(_omega/C);
        den3 = _epsilon[i]*k[i+1]/(_omega/C)+_epsilon[i+1]*k[i]/(_omega/C);
        c1 = complex<double>(cos(k[i+1]*d*um), sin(-k[i+1]*d*um));
        c2 = complex<double>(cos(_omega/par_v*d*um),sin(_omega/par_v*d*um));
        sf[i] = num1 * num2 / (den1 * den2 * den3) *c1 * c2;
        index = i % para.get_diel_n();
        d += l[index];
    }
}
//计算每层向后辐射的振幅
void source2(const double _omega, const double par_v, Para para, const double k[], const double _epsilon[], complex<double> sb[])
{
    int n = para.get_diel_n() * para.get_r_p_n();
    double l[para.get_diel_n()]; para.get_diel_l(l);
    int index;
    double d = 0.0;
    double temp = _omega * _omega / (C * C) - k[0] * k[0];
    double num1, num2; double den1, den2, den3;
    complex<double> c1, c2;

    for(int i = 0; i < n+1;i++)
    {
        num1 = par_v/C*(temp*C*C)/(_omega*_omega*_epsilon[i])*(_epsilon[i+1]-_epsilon[i]);
        num2 = 1-par_v*par_v/(C*C)*_epsilon[i]+par_v/C*k[i+1]/(_omega/C);
        den1 = 1-par_v*par_v/(C*C)*_epsilon[i]+temp*par_v*par_v/(_omega*_omega);
        den2 = 1+par_v/C*k[i+1]/(_omega/C);
        den3 = _epsilon[i]*k[i+1]/(_omega/C)+_epsilon[i+1]*k[i]/(_omega/C);
        c1 = complex<double>(cos(k[i]*d*um), sin(k[i]*d*um));
        c2 = complex<double>(cos(_omega/par_v*d*um),sin(_omega/par_v*d*um));
        sb[i] = num1 * num2 / (den1 * den2 * den3) *c1 * c2;
        index = i % para.get_diel_n();
        d += l[index];
    }

}


//计算前向出射的幅度
complex<double> amp_f(const int n, const complex<double> sf[], const complex<double> sb[], const complex<double> GRb[], const complex<double> _m_N[])
{
    complex<double> amp;
    amp = sf[n];

    for (int i = 0; i < n; i++)
    {
        amp += sf[n-1-i]*_m_N[i];
        amp += sb[n-i]*GRb[n-i-1]*_m_N[i];
    }

    return amp;
}


//计算后向出射的幅度
complex<double> amp_b(const int n, const complex<double> sf[], const complex<double> sb[], const complex<double> GRf[], const complex<double> _m_O[])
{
    complex<double> amp;
    amp = sb[0];

    for (int i = 0; i < n; i++)
    {
        amp += sb[i+1]*_m_O[i];
        amp += sf[i]*GRf[i+1]*_m_O[i];
    }

    return amp;
}

//粗略的积分函数
double my_int(double value[INT_N1][INT_N2], int row, int col, double dx1, double dx2)
{
    double integral = 0.0;
    for (int i = 0; i < row - 1; i++)
    {
        for (int j = 0; j < col-1; j++)
        {
            integral += (value[i][j] + value[i+1][j+1]) / 2.0 * dx1 * dx2;
        }
    }

    return integral; //尚存疑
}

//计算因子I_j
complex<double> facI(Para para, int j, const complex<double> GRf[], const complex<double> GRb[], const double k[])
{
    complex<double> tmp;
    double l[para.get_diel_n()]; para.get_diel_perm(l);
    int index = (j-2)%para.get_diel_n();
    complex<double> c(2*k[j-1]*l[index], 2*k[j-1]*l[index]);

    tmp = GRb[j-2] * GRf[j-1] * c;
    return tmp;
}
//计算各个域的因子I_j
void facI(Para para, const int n, const complex<double> GRf[], const complex<double> GRb[], const double k[], complex<double> fI[])
{
    for(int i = 0; i < n; i++)
    {
        fI[i] = facI(para, i+2, GRf, GRb, k);
    }
}

//计算域j内的前向电场分布
complex<double> Ezf(int j, const int n, const complex<double> sf[], const complex<double> sb[], const complex<double> GRf[], const complex<double> GRb[], const complex<double> _N[], const complex<double> _O[], const complex<double> fI[])
{
    complex<double> amp;
    complex<double> tmp[N]={0.0};
    complex<double> tmpf(0.0,0.0),tmpb(0.0,0.0);
    if (j == n+2)
    {
        relatedN(j, n, _N, tmp);
        return amp_f(n, sf, sb, GRb, tmp);
    } else if (j == 2)
    {
        tmpf = sf[j-2]; 
    } else if (j>2 & j<n+2)
    {
        relatedN(j, n, _N, tmp);
        tmpf = sf[j-2];
        for(int k = 1; k < j-1;j++)
        {
            tmpf += (sf[k-1] + sb[k]*GRb[k-1])*tmp[k-1];
        } 
    } else if (j == 1)
    {
        amp = 0.0;
        return amp;
    } else 
    {
        std::cout << "函数Ezf中参数j传递错误" << std::endl;
        return complex<double>(0.0);
    }
    
    for(int i = 0; i < N; i++)
    {
        tmp[i] = 0.0;
    }
    relatedO(j, n, _O, tmp);
    tmpb = sb[j-1];
    for(int k = j+1; k < n+2; k++)
    {
        tmpb += (sf[k-2]*GRf[k-1] + sb[k-1]) * tmp[k-j-1];
    }
    tmpf = tmpf * 1.0 / (1.0 - fI[j-2]);
    tmpb = tmpb * 1.0 / (1.0 - fI[j-2]) * fI[j-2] / GRf[j-1];
    amp = tmpf + tmpb;
    //std::cout << amp << std::endl;
    return amp;
}

complex<double> Ezb(int j, const int n, const complex<double> sf[], const complex<double> sb[], const complex<double> GRf[], const complex<double> GRb[], const complex<double> _N[], const complex<double> _O[], const complex<double> fI[])
{
    complex<double> amp;
    complex<double> tmp[N]={0.0};
    complex<double> tmpf(0.0, 0.0), tmpb(0.0, 0.0);

    if(j == 1)
    {
        relatedO(j, n, _O, tmp);
        return amp_b(n, sf, sb, GRf, tmp);
    } else if(j > 1 && j < n+1)
    {
        relatedO(j, n, _O, tmp);
        tmpb = sb[j-1];
        for(int k = j+1; k < n+2; k++)
        {
            tmpb += (sb[k-1] + sf[k-2] * GRf[k-1]) * tmp[k-j-1];
        }
    }else if(j == n+1)
    {
        tmpb = sb[j-1]; 
    }
    else if(j == n+2)
    {
        amp = 0.0;
        return amp;
    } else
    {
        std::cout << "EZb函数中参数j传递错误" <<std::endl;
        return 0.0;
    }

    for(int i = 0; i < N; i++)
    {
        tmp[i] = 0.0;
    } 
    relatedN(j, n, _N, tmp);
    tmpf = sf[j-2];
    for(int k = 1; k < j-1; k--)
    {
        tmpf += (sf[k-1] + sb[k] * GRb[k-1]) * tmp[k-1];
    }
    tmpf = tmpf * 1.0 / (1.0 - fI[j-2]) * fI[j-2] / GRb[j-2];
    tmpb = tmpb * 1.0 / (1.0 - fI[j-2]);
    amp = tmpf + tmpb;
    return amp;   
}
